Astronomy & Astrophysics manuscript no. code2DlD'astro-ph 


© ESO 2008 


February 5, 2008 





Simulating planet migration in globally evolving disks. 

A. Crida 1 , A. Morbidelli 1 and F. Masset 2 

1 Observatoire de la Cote d'Azur, B.P. 4229, 06304 Nice Cedex 4, FRANCE 
e-mail: crida@obs-azur.fr 

2 UMR AIM, DSM/DAPNIA/SAp, Orme des Merisiers, CE-Saclay, 91191 Gif/Yvette Cedex, FRANCE ; 
IA-UNAM, Apartado Postal 70-264, Ciudad Universitaria, Mexico City 04510, MEXICO 

Received : June 21, 2006 / Accepted : August 29, 2006 

ABSTRACT 

Context. Numerical simulations of planet-disk interactions are usually performed with hydro-codes that - because they consider only an 
annulus of the disk, over a 2D grid - can not take into account the global evolution of the disk. However, the latter governs planetary migration 
of type II, so that the accuracy of the planetary evolution can be questioned. 

Aims. To develop an algorithm that models the local planet-disk interactions together with the global viscous evolution of the disk. 

Methods. We surround the usual 2D grid with a ID grid ranging over the real extension of the disk. The ID and 2D grids are coupled at their 

common boundaries via ghost rings, paying particular attention to the fluxes at the interface, especially the flux of angular momentum carried 

by waves. The computation is done in the frame centered on the center of mass to ensure angular momentum conservation. 

Results. The global evolution of the disk and the local planet-disk interactions are both well described and the feedback of one on the other 

can be studied with this algorithm, for a negligible additional computing cost with respect to usual algorithms. 

Key words. Methods: numerical ; Accretion, accretion disks ; Solar system: formation 

locked in the middle of the gap, because both the outer and the 
inner part of the disk are repellent for it. Thus, the migration 
of th e planet has to go in parallel with t he migration of the 
gap jLin & PapaloizovJI 19861 1 Ward! 1200 3l) . which in turn has 
to follow the viscous evolution of the disk (characterized by 
a radial spreadi ng of the disk as the gas is a ccreted onto the 
central star (see lLvnden- Bell & Pringlel fT974i). This is usually 
referred to as type II migration. 

Although analytic theories have brought a great deal of 
understanding of the fundamentals of planet-disk interactions, 
it has become increasingly evident in the last years that nu- 
merical si mulations are an essential tool of investigation (see 
Pap aloizou et all2005l for the most recent review). Numerical 
simulations, however, are difficult and time-consuming, and 
they inevitably involve some simplifications, which might turn 
out to produce artifacts. In this paper we are concerned with the 
simulation of type II migration. 

As we have seen above, a proper simulation of type II 
migration requires not only the correct calculation of planet- 
disk interactions - which are essentially local - but also of the 
global evolution of the disk. Unfortunately, planet migration is 
typically simulated with hydro-dynamical codes using a 2 di- 
mensional polar grid which, for numerical reasons, is truncated 
at an inner and an outer radius. This enables to describe well the 



1. Introduction 

Planetary formation occurs in disks of gas and dust around 
protostars. In particular, giant planets - whose mass is mostly 
made of hydrogen and helium - must have formed before the 
dissipation of the gas disk. Consequently, they must have ex- 
erted tidal forces on the gas and their orbits must have evolved 
in response to the gas. 

In particular, the presence of a planet on a circular orbit 
lea ds to the formation of a s piral density wake in the disk . 
iGoldreich & Tremaind dl98Ct) and llin & Papaloizou! lll979h 
showed that, through the wake, the planet gives angular mo- 
mentum to the part of the disk exterior to its orbit, while it 
takes a ngular mom entum from the inner part. For small mass 
planets, War dl lll997l) showed that the net result is a loss of an- 
gular momentum for the planet, which makes its orbit decay on 
a short timescale. This is usually referred to as type I migration. 

As the planet gives angular momentum to the outer part 
of the disk, it repels it outward; symmetrically it repels the 
inner part inward. If the planet is massive enough (the threshold 
mass depending o n the disk's viscosity and scale height ; see 
ICrida et all J2006ll and references therein) a clear gap opens 
around the planet's orbit, effectively splitting the disk into an 
internal and an external part. In this situation, the planet gets 
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local interaction of the planet with the disk but not the global 
evolution of the disk. 

Indeed, the boundary conditions at the extreme rings of the 
polar grid can not take into account what happens in the whole 
disk, outside of the grid. For instance, the use of open bound- 
ary conditions allows the gas to leave the region covered by 
the grid ; this makes the studied disk annulus behave as if it 
were surrounded by vacuum, so that it empties very rapidly, 
which is obviously not realistic. In the opposite extreme case, 
boundary conditions which impose that the mean density on 
the extreme rings remains constant with time, disable the ac- 
cretion and spreading of the gas. Some other prescriptions be- 
tween these two extremes may be used, (for instance allowing 
inflow, or imposing an outflow given by an analytical model, or 
setting the flow on the last ring equal to the viscous flow mea- 
sured in the neighboring rings, etc.). However, it is very diffi- 
cult to adapt these prescriptions to the changing behavior of the 
disk with time, in particular when the disk undergoes perturba- 
tion from the planets. In any case these prescriptions are rather 
arbitrary and, thus they may introduce possible artifacts in the 
planetary evolution. 

In this paper we present a novel idea for the correct cal- 
culation of the global evolution of the disk. It consists in sur- 
rounding the 2D polar grid with a ID radial grid. The ID grid 
extends from the real inner edge of the gas disk (e.g. the X- 
wind truncation radius at a few tenths of AU) to the real outer 
edge (e.g. the photo-dissociation radius at hundreds of AU). 
This ID grid has open boundaries at the inner and outer edges, 
and exchanges information with the 2D grid for the definition 
of realistic, time-dependent boundary conditions of the latter. 
Our algorithm for the interfacing between the ID and 2D grids 
is driven by the requirement that the angular momentum of the 
global system (the disk in the 2D section, plus the disk in the 
ID section plus the planet-star system) is conserved. We will 
describe it in detail in Sect. [5] First, however, we need to revisit 
the algorithm used for modeling the gas evolution and planet- 
disk interactions on the 2D grid, to ensure that it also conserves 
angular momentum, which is often not the case in standard im- 
plementations. We discuss this issue in Sect. |2] In Sect.|4] the 
results of our new algorithm will be discussed, in terms of CPU 
time and robustness with respect to the positioning of the inter- 
faces between the two grids. Finally in Sect.|5j we will describe 
some interesting astrophysical applications of this new hydro- 
dynamical code. 

2. Conservation of angular momentum in 2 
dimension hydrodynamical algorithms 

Consider a simulation of the planet-disk interactions, with the 
disk represented on a 2D polar grid with open boundary condi- 
tions. A necessary requirement for the simulation to be correct 
is that the sum of the angular momenta of the disk, of the star- 
planers) system and of the gas outflowing through the bound- 
aries, remains constant over time. 

We have tested if this is the case, using the code FARGO 
(lMassel2000l) . 

Our simulation accounts for a Jupiter mass planet on an ini- 
tially circular orbit. Here and in the rest of the paper we adopt 
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Fig. 1. Relative variation of mass (top panel) and angular mo- 
mentum (middle panel) of the system {gas + Jupiter + outflow) 
with time. The absolute variation of angular momentum of the 
system is compared to that of the planet on the bottom panel (in 
our units, in which the initial angular momentum of the planet 
is 1(T 3 ). 



the following units : the solar mass, the initial semi-major axis 
of the planet and its orbital frequency, so that the gravitational 
constant G = 1 and an orbit lasts 2n time units at r = I. The 
grid used to represent the disk extends from r = 0.25 to r = 3 
with open boundaries. It is equally divided in N r = 165 elemen- 
tary rings, and N s = 320 sectors. The disk aspect ratio (H/r) 
is set uniformly to 5%. It is assumed to be constant in time, 
hence the disk is locally isothermal. The equation of state used 
in FARGO is : P = c, 2 2, with c s - HQ. as usual, where P is 
the pressure, 2 the density, c s the sound speed, H the disk scale 
height, and Q the local angular velocity. The gas kinematic vis- 
cosity is v = 10 • in our normalized units (which corresponds 
t o the viscosity at the location of the planet for a = 1 .25 10~ 3 in 
a 
S 



Shakura & Sunvae vi lli 9731) prescription). Its mean density is 



3.10 4 , which is a bit less than the Hayashi Minimal Mass 
Solar Nebula at Jupiter dHavashill98ll) ; the initial density pro- 
file can be seen in Fig.[3]and corresponds to a disk that evolved 
for some time under the effect of its own viscosity ; it can be 
very well approximated by E(r) = 0.000306 exp(-r 2 /52.8). 

The evolution of the total mass and angular momentum of 
the whole system (including the outflow, e.g. the cumulated 
mass and momentum advected outside of the region spanned 
by the 2D grid by the gas outflowing through the boundaries) 
is presented in Fig. [2 As one sees, while the total mass is con- 
served at the level of numerical errors (top panel), the conser- 
vation of the angular momentum is quite poor (middle panel). 

The bottom panel shows that the gain of angular momen- 
tum of the whole system (plain line), which is the error of the 
simulation, amounts to 10% of what the planet exchanges with 
the gas in its migration process (dashed line). Thus we can ex- 
pect that the migration of the planet measured in the simulation 
is correct only at the 10% level. 

This shows the need to improve the algorithm in order to 
achieve a much better conservation of angular momentum. We 
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describe a procedure that considerably improves angular mo- 
mentum conservation below. 

2.1. Choice of reference frame 

Need of a frame centered on the center of mass. Most, if not 
all, the published numerical simulations of disk-planet interac- 
tion that we are aware of used a non-inertial f rame ce ntered 
on the primary JCrida et alJl2006l iKleyl 1 1 999t iMassell 1 2002 ; 
Masset et al1l2006l jNelson et al.ll200(i iNelson & Bend £003 ; 
Varni ereetal.ll2004l l2005). This is done for practical reasons, 
and because it is thought to be better adapted to describe the 
motion of the inner part of the disk. This reference frame is not 
inertial, so that indirect forces are taken into account for the 
simulation to be realistic. It is well known that the presence of 
these forces destroys the conservation of the angular momen- 
tum measured in the non-inertial frame. In principle however, 
one can make at any time a change of coordinates to compute 
position and velocities in the inertial frame (centered on the 
center of mass), and the angular momentum computed from 
these coordinates should remain conserved. This is for instance 
what happens in N-body codes, when the simulation is done in 
heliocentric coordinates (except from truncation and round-off 
errors). This is not the case here. The transport algorithm for 
the gas - in which mass, angular momentum, linear momentum 
etc. advect from a cell to its neighbors (sub-step 5 in Sect. 12. 2> 
- imposes the conservation of each of these quantities. This is 
correct for the mass, but not for the momenta, because they are 
not supposed to be conserved in the adopted frame. As a con- 
sequence of this imposed, unphysical conservation, the conser- 
vation of the momenta in the inertial frame is corrupted. The 
solution out of this problem is the use of the frame centered on 
the center of mass throughout the algorithm. 

Implementation. This is not trivial to implement. A first pos- 
sibility is to suppress all indirect forces in the algorithm and, 
whenever the position and the velocity of the star are needed, to 
compute them imposing that the position of the center of mass 
of the whole system (star, planets and disk) is at (0,0). This is, 
for instance, what is usually done in N-body barycentric simu- 
lations. In this case, however, the situation is more complicated. 
While in N-body codes all interactions are treated simultane- 
ously, planet-disk interaction simulators involve a purely grav- 
itational part and an hydrodynamical part, and treat the plan- 
etary system and the disk separately (see !2.2l for the sequence 
of integration sub-steps). When the disk is advected, the star is 
moved to assure the position of the center of mass, but not the 
planet. Thus, at the next step, the planet sees a star in a differ- 
ent location in phase-space with respect to the end of the previ- 
ous step. This has dramatic effect for the orbital stability of the 
planet and the overall conservation properties. A second, more 
advantageous possibility is to consider the star exactly like a 
planet during all the stages of the computation. This ensures 
a symmetry in the treatment of the star and of the planet(s). 
When the gas is advected, the planetary system (now including 
the star) is translated in phase-space to ensure that the center 
of mass of the whole system remains motionless at r — 0. This 



little correction does not perturb the relative planetary motion. 
It corrupts, in principle, the conservation of the total angular 
momentum, but the errors introduced have negligible conse- 
quences, as we will test and discuss below. 

Caveats. Because the grid is centered on the center of mass, a 
good numerical accuracy can be achieved only if the motion of 
the gas is approximated by a rotation around the center of mass. 
This is the case as long as the central object is close to the cen- 
ter of mass, namely if the total mass of the planets is negligible 
with respect to the one of the star (which is the common case) 
or - in case of a large stellar companion - if one considers a 
distant circumbinary disk. The choice of a frame centered on 
the center of mass is obviously not adapted for studying a disk 
around a star with a massive companion or a circumplanetary 
disk. 

In the case of an axisymmetric problem (accretion disk with 
no planet), the choice of such a frame may lead to numerical 
instabilities because small deviations from axial symmetry in 
the disk due to numerical errors cause a shift of the star rela- 
tive to the center of mass, which in turn can enhance the disk's 
asymmetries. Thus, if the planetary system is made of the sole 
star, we recommend to keep it fixed at the origin of the frame. 

Finally, we remark that in many applications, the equations 
of motion are implemented in a rotating frame in order to keep 
the planet at a constant position. This can still be done in a 
frame centered on the center of mass. 

2.2. Sub-steps sequence 

As we anticipated above, the integration algorithm treats sepa- 
rately the pure gravitational part and the hydrodynamical part. 
To ensure a good preservation of the conserved quantities, it 
is necessary to respect as much as possible the action-reaction 
principle in the planet-disk interactions. This requires that the 
sequence of integration sub-steps is taken in a specific order. 

Here is the sequence that we adopt : 

1 . The gravitational potential of the planets and star is com- 
puted and stored. 

2. The velocities of the planets and star are updated using the 
gravitational influence of the disk. 

3. The velocities of the gas are changed in each cell according 
to the non advective part of the Navier-Stokes equations 
(the external forces dues to the gradients of pressure and 
gravity field and the viscous stress) ; the gravitational po- 
tential computed at sub-step 1 is used. 

4. The planet-star system is advanced under its own gravita- 
tional interactions using a N-body algorithm (specifically a 
5th order Runge-Kutta integrator). 

5. The advective part of the Navier-Stokes equations is per- 
formed: using the disk velocities computed at sub-step 3, 
the mass, angular and radial momentum are advected from 
each cell to the neighboring cells. The new velocities of the 
gas are then computed in each cell from the new momenta 
and masses. 
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6. The conservation of the center of mass is ensured by trans- 
lating in space and velocity the center of mass of the 
planets-star system. 

As one sees, in the above algorithm, every time that an action 
is made on the planetary system, the equivalent action is im- 
mediately applied on the disk. Moreover, there is no difference 
in treatment between planets and star, which is essential for 
the reasons discussed above. However, particular care has to 
be paid in the implementation of the algorithm to ensure the 
effective symmetry of these actions ; this is discussed in next 
subsection. 

2.3. Action-reaction symmetry 

The action-reaction principle has to be perfectly fulfilled in the 
computation of sub-steps 2 and 3. Notice that it is not the case, 
for instance, if one computes the action of a planet on the disk 
from the gradient of its potential on each cell, and the action 
of the gas on the planet from the sum of the elementary grav- 
itational forces that it feels from each cell. In principle, both 
calculations should be equivalent, giving total forces of the 
same strength and opposite directions. However, if the gradient 
is computed by finite differences (typical in grid calculations), 
this introduces a difference with respect to the force exerted by 
the corresponding cell on the planet. 

In order to impose action-reaction symmetry, we proceed as 
follows. After sub-step 1, we measure the total change of an- 
gular momentum that the gas will have at sub-step 3 due to the 
potential of each planet. Then, in sub-step 2, we compute the 
change of azimuthal velocity of each planet according to the 
change of angular momentum that it causes to the disk. This 
ensures angular momentum conservation between sub-step 2 
and sub-step 3. Similarly, in sub-step 2 the planet's radial ve- 
locity change is computed by measuring the total change of 
linear momentum in this direction caused to the gas disk by the 
same planet. 

2.4. Results and discussion 

All these precautions allow us to improve greatly the angular 
momentum conservation of the whole system. The same sim- 
ulation of Fig.^is computed with our modified algorithm. We 
denote by 6H the variation of angular momentum, in the iner- 
tial frame centered on the center of mass, of the whole system 
(the sum of the momenta of the gas in the 2D grid, of the planet- 
star system, and of the outflow which, we remind, is the cumu- 
lated momentum carried by the gas that flew out of the grid). 
It should be zero if the code were perfectly conservative. Thus 
5H is a measure of the error of the integration scheme. It is 
compared with three relevant quantities in Pig. El The bottom 
curve (solid line) corresponds to the logarithm of the relative 
variation of angular momentum of the system: SH/Hq, where 
Ho is the initial total angular momentum (it corresponds to the 
middle panel of Fig. [0 ; on a long timescale (16000 time units, 
> 2500 orbits at r — 1), this normalized error does not exceed 
10~ 5 5 . The middle curve (long dashed line) corresponds to the 
logarithm of SH/H p , where H p is the angular momentum of the 











- 


log (|d H / H_p | ) 






ill / 


log ( | d H / H_0 | ) 


I 



2000 4000 6000 8000 10000 12000 14000 16000 



Fig. 2. Relative angular momentum variation as a function of 
time. The variation of angular momentum of the system { gas 
+ planet + star + outflow }, 5H, is compared to its initial value 
//o (bottom curve, plain line), to the angular momentum of the 
planet H p (middle curve, long-dashed line), and to the variation 
of the latter 5H p (top curve, short-dashed line). 

planet : the error in total angular momentum is about 4.5 orders 
of magnitude smaller than the angular momentum of the planet. 
Thus, it should not affect its migration. Indeed, the top curve 
(short dashed line) shows the logarithm of the ratio between 
5H and 6H p , where 6H p is the variation of angular momentum 
of the planet : this ratio is smaller than 10" 3 . 

Although highly satisfactory for our purposes, the conser- 
vation is admittedly not perfect, in particular when compared to 
the mass conservation (which is in the new simulation as good 
as in the example of Fig.^. This is due to two sources of error. 

The most obvious one is that the algorithm used for the 
advancement of the planetary system, which is in our case 
a Runge-Kutta algorithm, is not symplectic. However, we 
checked that in our simulation, the amount of angular momen- 
tum artificially introduced in the system in this way is 3 orders 
of magnitude smaller than the total global error. 

The second source of error is in sub-step 6, due to the trans- 
lation of the center of mass of the planet-star system. This 
translation occurs for two reasons, (i) It is required to compen- 
sate the center of mass motion due to errors introduced by the 
discretization of the grid and to the second order (in time) ad- 
vection of the disk. Thus, the amplitude of this effect decreases 
with increasing disk resolution and decreasing time-step, (ii) 
Gas continuously leaves the grid through its inner or outer 
boundary. This gas is deleted from the simulation, while its 
total mass and momentum are recorded as outflow. The prob- 
lem is that the outflowing gas is in general not axisymmetric. 
Therefore, its elimination from the simulation causes a shift (in 
position and velocity) of the center of mass of the whole sys- 
tem. The shift is artificial, in the sense that it would not exist if 
the gas were modeled with an infinite 2D grid, and it requires, 
for compensation, the translation of the planetary system. This 
issue is a conceptual one. Thus this source of error cannot be re- 
duced by tuning the grid resolution or time-step. We think that 
there is no way around it : it is the price to pay to work with 2D 
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grids that are not as extended as the physical disk. However, as 
we have seen in Fig. [2 this noneliminable error is so small that 
for all practical purposes we can safely ignore it. 

3. Coupling of the 2D grid with a 1 D grid 

As we discussed in the introduction, for a correct simulation 
of type II migration the conservation of angular momentum is 
a necessary, but not sufficient condition. It is also necessary 
that the global evolution of the disk is correctly reproduced, 
which is not the case if only a portion of the disk (the annulus 
represented by the 2D grid) is considered. 

For this reason, we consider the whole disk, from its real 
inner to outer radius, c orrespondi ng respectively to the 'X' ra- 
dius at the inner edge JShull20Q2|). and to a photo-eyaporatiq n 
radius at the outer edge (see Hollenbach & Adams 2004a b). 
Because the representation of the whole disk with a 2D grid 
would be computationally prohibitive, outside of the radial dis- 
tance range covered by the usual 2D grid we study the evolu- 
tion of the extended disk in one dimension, with the aid of a 
ID grid representation. In practise, the ID grid is made of two 
parts : one extending from the real inner edge of the disk to the 
inner edge of the 2D grid, and the other one extending from the 
outer edge of the 2D grid to the real outer edge of the disk. 

In this section, we first present the equations for the ID 
disk. Then we detail how we couple the ID and 2D portions 
of the disk in a way that conserves mass and angular momen- 
tum, and that ensures a proper simulation of the disk's global 
evolution. 

3.1. 1D equations 

The Navier-Stokes equation and the advection of the gas are 
computed in the ID grid exactly like in the 2D grid, but with 
N s = 1, so that each cell can be considered as an elementary 
ring. Consequently, there are no azimuthal components for the 
various gradients, neither an advective transport phase in the 
azimuthal direction. The gravitational potential due to the plan- 
ets and the star is computed for the ID grid as a function of r 
only : the whole planetary system is considered as a single cen- 
tral mass, the mass of which is the sum of all the bodies inside 
the considered radius. Denoting by fi the Heaviside function 
equal to 1 for positive arguments and to for negative ones, the 
gravitational potential felt by a cell of the ID grid located at a 
radius r from the center of mass reads : 

, n v-i GM,, , 
* M = 2j -T^ir-re) 

p 

where the index p goes through the whole planetary system 
(including the star). This means that the potential felt is equal 
to —G/r multiplied by the mass of all the celestial bodies that 
lie inside the ring. Thus the position of the star with respect to 
the center of mass has no influence ; it would be the same if the 
ID rings were centered on the star. 

We assume that the gas disk in the ID grid is axisymmetric, 
so that in principle there is no torque in the interaction between 
the planets and the ID grid. For more realism, however, one 
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Fig. 3. Sketch of the coupling via ghost rings at the inner inter- 
face between the ID grid and the 2D grid. 



can assume that each ID ring fe els a gr avitational torque due to 
ea ch planet, using the formu la o flGoldreich & Tremainel(ll980l) 
or lLin & PapaloizouN 19791) : 

5T g (r) * 0.4 ^jZr/O/r" 1 {^j (2nrY.5 r ) (1) 

where A = r - r p , Q p is the angular velocity of the planet, and 
6r is the width of the ring. Of course, for angular momentum 
conservation, the opposite of this torque exerted on the ID grid 
has to be exerted on the planet. 

Because of the assumption of axial symmetry, what leaves 
the ID grid through its boundaries does not impose any re- 
adjustment of the star-planet(s) center of mass. Thus, the con- 
ceptual problem discussed at the end o fl2.4l is not relevant here. 
Consequently, the evolution on the ID grid always shows a per- 
fect angular momentum conservation, at the level of numerical 
errors. 

Despite of this perfect conservation, however, for the re- 
sults to be correct it is necessary that the assumption of axial 
symmetry is a good approximation for the real disk. Moreover, 
because Eq. Q is an approximation, it is necessary that the 
planetary torque is small for the error to be small in absolute 
value. Both requirements are fulfilled if the interfaces between 
the ID and 2D portions of the disk are placed sufficiently far 
from the planet(s). We will detail on this issue in Sect. |4] with 
quantitative tests. 

3.2. Ghost rings 

The ID grid and the 2D grid communicate with each other 
through a system of ghost rings. This technique is inspired from 
that used for the parallelization of a hydro-code on a distributed 
memory architecture. We describe it here for the inner interface 
between the ID grid and the 2D grid, that is the interface at the 
inner edge of the 2D grid ; the ID grid is inside, extending from 
the real inner edge of the disk (the truncation radius or the sur- 
face of the star) to the 2D grid ; the 2D grid lies outside this 
interface, as sketched in Fig. [5] 

Outside the interface, a number Nc of ID rings are added, 
which superpose to the first Nc rings of the 2D grid ; they are 
the ghosts for the inner ID grid (see Fig. [5}. At the beginning 
of every time-step, the density and the velocity in the ID ghost 
rings are set as follows. The density of gas (S), radial momen- 
tum (Evy), and angular momentum (Ervg) in the ghost ring are 
set to the azimuthal average of the respective quantities in the 
corresponding ring of the 2D grid. Denoting with the super- 
script 'ID' the quantities in ID rings and by h the specific an- 
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gular momentum (h = rvg), this leads to the following equa- 
tions : 



S 1D v,. 1D 
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ID 



i r 2 " 

Zrd6>, 

2?rr Jo 



- r- f : 



Sv r rd9 . 



2nr 



2k 



rd6 . 



(2) 
(3) 
(4) 



As E' u is given by Eq. 0, v r and vg 1D are easily computed 
from Eq. Q and @. 

Then, during the time step, the computation of all stages is 
performed in the ID grid as well as in the ghost rings. The ID 
grid 'feels' the ghosts via the Navier-Stokes equations. As the 
ID ghost rings have been filled with quantities inherited from 
the 2D grid, the ID grid behaves as if it felt the 2D grid outside 
the interface. 

Symmetrically, inside the interface, a Nc number of 2D 
rings are added : the 2D grid ghosts. They are treated as normal 
2D rings, except that at the beginning of every time step, the 
azimuthal means of surface density, of radial momentum and 
angular momentum density are set equal to the respective val- 
ues measured in the corresponding ID rings, as sketched with 
arrows in Fig. [3] So, the 2D grid effectively 'feels' the ID grid 
inside the interface. Filling the 2D ghosts is a bit more elabo- 
rated than in the ID ghost rings case. Indeed, for a given 2D 
ghost ring, one has to proceed in four steps : (i) store for each 
cell the surface density, the radial momentum, and the angular 
momentum ; (ii) compute the azimuthal means of these quan- 
tities ; (iii) subtract in each cell, for each of the three consid- 
ered quantities, the difference between its azimuthal mean and 
its value in the corresponding ID ring; (iv) find the velocity 
in each cell by dividing the new momenta by the new surface 
density. 

The minimal number of ghost rings needed, Nc, depends 
on the numerical scheme : it is the number of rings that causally 
affects a given ring during a time-step, or equivalently the num- 
ber of rings on which the information contained in a ring prop- 
agates during a time-step ; it is often called the kernel. For in- 
stance in FARGO, N c = 5. 

This way of coupling the ID and 2D grids ensures a smooth 
transition for each of the quantities computed in the code. 
In particular, if there is no planetary perturbation, the global 
disk (2D + ID parts) behav es exactly as predicted by the 
iLvnden-Bell & Pring ld i 19741) equations. This can be seen in 
Fig.|4] which corresponds to Fig. 4 in the aforementioned pa- 
per : it shows the evolution of the density distribution at four 
time s, for £ >= p(r) _ 6(r- 1), with 6 the Dirac distribution. Like 
m lEvn den- Bell &PringleNl 974 . T * denotes the viscous time : 
r* = 6vt in our units. The little discrepancy between the theo- 
retical profile (thin) and the numerical one (bold) at T* = 0.004 
most likely comes from the fact that the initial distribution is 
not exactly a 6 function; this little difference vanishes with 
time. 




0.8 1.2 1.6 

distance to central star 

Fig. 4. Viscous spreading of a ring. The theoretical 
curves (obtai ned by numerical solution of Eq. (14) of 



Lvnden-Bell & Pringy ill 9741) with the boundary condition 
2(0.02) = 0) are thin, while the numerical ones obtained with 
our code are bold. The vertical dashed lines show the position 
of the interfaces, where the solution provided by our code 
remains perfectly smooth. 

3.3. Fluxes at the interface 

The fluxes of mass and angular momentum through the inter- 
faces have to be perfectly conservative. What leaves the 2D grid 
has to enter the surrounding ID grid, and vice-versa. 

3.3.1. Mass flux 

The definition of v r 1D in Eq. makes the mass flux computed 
in ID be the same as the one computed in 2D. Indeed their 
theoretical expressions write : 



2D 



ID 



2k 



"Ev r rd6 



o 

= 2nfL 



ID ID 

1 r 



One can see that Eq.|3]is equivalent to F m 1D = F m 2D . However, 
in a numerical scheme with finite time steps and discrete grid, 
the flux is most often not computed so straightforwardly. Care 
has to be taken that the numerical expressions of the flux coin- 
cide at the interface between the ID and the 2D grid. In stag- 
gered mesh codes, the density and the velocity are not defined 
at the same locations in the cells of the grid ; typically, the den- 
sity is given at the center, while the radial component of the 
velocity is given on the middle point of the inner edge and the 
tangential componen t on the middle p oint of the left edge (see 
for instance Fig. 3 of S tone & Normanll993t) . Thus, the fluxes 
are most often computed as : 

r ID _ ylD* ID 

1 m computed — ' - £j v r 

^m 2D computed = ^J^" [j]v r [j] r69 

j 

where the index j goes through all the cells of the ring, 56 is the 
angular size of a cell, r is the radius of the interface, and E* and 
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£ 1D * are the values o f the densities £ and £ 1 D at the interface 
after half a time-step JStone & Normanl 19931 paragraph 4.4). 

To ensure the mass conservation, we are led to solve the 
following equation instead of Q for the ring next to the inter- 
face : 



7 ID 
m computed 



= F, 



2D 



computed 



(5) 



This equation is only slightly different from Eq. Q. It is 
an implicit equation in V r 1D , because £ 1D * depends on v r lD . 
Nevertheless, it can be solved by iteration. 

3.3.2. Angular momentum flux 



Denoting with a prime sign the non axisymmetric components 
of the quantities (h! = h — h lD , v' r — v r - v r 1D ), using Eq. ©-(O 
the angular momentum fluxes in ID and in 2D write: 



F h 1D = 2nrl} D h™v r m 



£/zv r 





r&O = F h m + 



Zh'v' r 

Jo 



rd6 



Independently of the numerical scheme and of the modifica- 
tion of v r 1D seen before, it appears that it is impossible for 
the ID flux to equal the 2D flux. While F/, 1D corresponds to 
the angular momentum carried by the gas flowing through the 
interface, the term F' h = f Yh'V r rd6 corresponds to a flux 
of angular momentum that is not due to advection by the ax- 
isymmetric flow ; it comes from the azimuthal perturbations 
of v r and Vg, that represents the angular momentum carried 
by a wave. The wave-carried propagation of angular momen- 
tum is well-kn own in planet-disk interactions (see for instance 
Appendix C oflCrida et alJ2006llGoldreich & Nicholson! 19891 
iTakeuchi et al J 19961) . but cannot be accounted for by a ID grid. 
This is a conceptual problem. 

In fact, we face a degree of freedom problem. With only 3 
free variables (£ 1D , v r lD , and v# 1D ), 4 quantities have to be set 
in the ID ghosts : the mass and the angular momentum that is 
in each ring, the mass and angular momentum fluxes. The first 
two are set with the prescriptions Q and ®, which determine 
the values of £ 1D and vg lD . The third variable, v r lD is naturally 
set by Eq. 0/Eq. (|5j- It appears here that it is impossible to 
obtain the correct flux of angular momentum in the ID grid 
ghosts. 

A possible solution is that, during the advection phase in 
the ID part of the algorithm, one imposes that the angular mo- 
mentum flowing through the interface corresponds to the flux 
computed in the 2D grid, f a 2D computed- However, as the flux of 
angular momentum carried by a non axisymmetric wave cannot 
be represented with ID Navier-Stokes equations, F', would be 
deposited abruptly in the first ring of the ID grid. This would 
lead to the formation of a spurious gap at the interface, and 
possibly lead to a numerical instability. 

It has been shown that pressure supported waves travel 
through the disk and deposit their angular mome ntum 
smoothly, as they g et damped viscously jTakeuchi et ai1 ll996> 
or through shocks ( Goodman & Rafi kovKOOll) . Thus, a better 
solution for our problem is to model this wave propagation 
through the ID grid. We first perform advection as usual in 



each grid and evaluate on the interface F' h = Fh - F/, (or, 
better, 5F h = F A 2D computed - F/, "Computed ~ F' h , which might be 
slightly different from the former according to the numerical 
scheme). Then, over a time-step 5t, we spread in the ID grid 
the amount of angular momentum 5F 6t. A prescription for 
the dep osition of the flux can be found in lGoodman & Rafikovl 
) as a function of the distance from the planet. This could 
be used when there is only one planet, but not if there are sev- 
eral planets, because it is difficult to know which fraction of 
6Fh is due to each planet. Thus, we adopt an exponential func- 
tion, because it is scale free so that, once 6F/, is known at the 
interface, its deposition does not depend on the position of the 
planet(s) and is a function of the distance from the interface 
only. In practice, we assume that, over a time-step 6t, the an- 
gular momentum deposited by the wave in a ring of width dr 
located at a distance d from the interface between the ID and 
2D grids is : 



6h = 5Fh 6t A exp( — ) dr . 

A 



(6) 



where A is the damping length-scale of the flux. This assump- 
tion will be supported in Sect. 14.11 where the evolution of the 
disk in the ID grid is found to be insensitive to the position of 
the interface with respect to the planet. 

We assume that A = 0.5 in natural units, which is the 
value obtained by fitting with an expo nential the decay o f the 
wave-carried flux in appendix C of ICrida et alJ J2006I) . For 
simplicity, we assume that A does not depend on the disk's 
viscosity and scale-height ( an assumption partially justified in 
iGoodman & Rafikovl feOOlk 

The deposit of the quantity dh of angular momentum in 
a ring is simulated by applying a suitable torque, namely by 
adding to v l g D the quantity 6v l g D = dh/r, where r is the radial 
distance of the ring from the star. 

Notice that the integral of l|6} from d = to infinity is equal 
to the total angular momentum 6Ff, 5t carried by the wave at 
the grid interface. However, the ID disk is not infinite in radial 
extent. Thus, a fraction of the angular momentum carried by 
the wave will not be deposited in the disk, but will outflow 
from the system. This outflowing momentum, as well as the 
angular momentum and mass advected through the inner and 
outer radius of the ID grid are recorded, in analogy with what 
was done in Sect.[2]for the 2D grid alone. 



3.3.3. Mean viscous torque 

The coupling of the ID and the 2D grid described in the two 
subsections above ensures a smooth, conservative evolution of 
the disk. The gas is free to accrete and spread from the star 
to its physical outer edge, through the interfaces between the 
grids. However, at the interfaces, the shear gives an azimuthal 
viscous stress. This appears as a torque exerted on a grid by its 
ghosts (and reciprocally). The torque exerted by the part of the 
disk inside a given radius ro on the part of the disk outside ro 
reads, from the Navier-Stokes equations : 



t v (ro) 



pin 

= r I : 
Jo 



T rB r dd , 
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where T = „ rr „ r is the stress tensor, T = 
\i0r lee ) 

TLv(p - (|Vv)/) , with D the strain tensor and 7 the identity 

matrix. Thus, T r6 = T Br = £v(±§£ + r^M). In a staggered 
mesh scheme, T r g is defined at the inner edge of every cell. 

For angular momentum conservation, the torque felt by 
the ID grid from its ghosts must equal the one exerted by 
the 2D grid on its ghosts. This requires T^P (/"interface) = 
3^ Zj TYeUKnnterface), where the index j goes through the N s 
cells of the 2D ring next to the interface. However, given that 
the expression of T r g implies the product of velocity gradients 
by the density, its average is not equal to the product of the av- 
erages, and the required equality is not necessarily true. Thus, 
one has to replace the value of TIP at the interface by the aver- 
age of T r g on the interface in the 2D grid. 

3.4. Results and discussion 

We consider the final state of the simulation of Sect.|2]: a Jupiter 
mass planet initially at r p = 1 evolves in a gas disk represented 
by a 2D grid extending from r = 0.25 to 3 for 16000 time 
units (» 2500 orbits). The final density profile of the gas disk 
is shown in Fig.|5]as a dashed line. One can see that the planet 
has opened a gap and migrated inward (the planet is located 
in the middle of the gap). Moreover, the surface density of the 
gas over the considered range is strongly reduced relative to its 
initial value (dot-dashed line in Fig. |5}. We then compute an- 
other simulation with the same planetary system and the same 
2D grid, but introducing a ID grid extending from 0.1 17 to 20 
length units (approximately from 0.5 to 100 AU, assuming the 
unit length equal to 5 AU). The final profile of the gas distribu- 
tion obtained in this new simulation is shown as a solid line in 
Fig. |5] Its bold part corresponds to the 2D portion of the disk. 

We remark two important aspects in Fig. [5] First, the radial 
profile of our new solution is clearly smooth, which indicates 
that there are no artifacts in the passage of information from the 
2D grid to the ID grid and vice-versa. The small kink visible at 
the inner boundary of the ID grid (r ~ 0.15) is due to the imple- 
mentation of the open boundary condition (see appendix |A) ; 
this artifact also appears in Figs.|4][7]and[8] The change of sign 
of the second order derivative of the radial profile near the outer 
interface (r ~ 3) is a real feature at the considered time-step, 
as we will discuss later (see Fig.|9]and the related discussion). 
Second, despite the disk in the new simulation has also signif- 
icantly evolved relative to the initial profile, the new surface 
density of the disk is very different from that obtained in the 
simple 2D simulation (compare the solid and dashed curves). 
In particular, in the new simulation the outer part of the disk 
has not significantly evolved, which is due to the long viscous 
time there, as the latter scales with r 2 /v. Conversely, in the old 
simulation, because of the open boundary condition, we ob- 
served a significant disk erosion. This shows the importance of 
the boundary conditions for the global evolution of the system. 

The global conservation in the system (namely the evolu- 
tion of the mass and angular momentum of the planetary sys- 
tem, plus those of the gas in the two grids, plus those advected 
through the inner and outer radius of the ID grid and recorded 



0.00035 




distance to central star 

Fig. 5. Gas surface density profile after 16000 time units ~ 
2500 orbits of the Jupiter mass planet initially placed at r p = 1 . 
The solid line corresponds to a 2D (bold) and a ID (thin) grid 
coupled, while the dashed line stands for a 2D grid alone. The 
dot-dashed line shows the initial profile. 



1 D and 2D grid coupled 2D grid alone 
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time time 



Fig. 6. Relative variation of mass and angular momentum in the 
case of a 2D grid with open boundaries and a 2D grid coupled 
with a ID grid. 

as outflow) is presented in Fig.|6j and compared to the conser- 
vation obtained in the code using the 2D grid only. One can see 
that the use of the extended ID disk does not change the excel- 
lent result obtained in Sect. |2 Thus, the coupling between the 
two grids that we described above is correct in terms of con- 
servation properties. Actually, the error in angular momentum 
is slightly larger in the case with the ID grid because the den- 
sity at the inner interface of the 2D grid does not vanish, as can 
be seen in Fig. [5]; so, the major source of error (the non ax- 
isymmetric outflow from the 2D grid, discussed at the end of 
Sect. 12.41 does not disappear, while it does for the simulation 
using the 2D grid alone. 

4. Performance of the code 

In this section, we check the accuracy of the results, as a func- 
tion of the size of the 2D grid, and we discuss the CPU cost 
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Fig. 7. Inner disk profiles at the end of simulations with differ- 
ent locations of the inner interface between the ID and the 2D 
grid. The bold part of the profiles corresponds to the 2D grid. 

of this hybrid scheme. We show that the coupling of a ID grid 
increases by far the realism of the results of the hydro-code, for 
a low additional computation cost. 

4.1. Position of the interfaces with respect to the 
planet. 

As already mentioned it in sec. l3.ll the interface between the 
2D and the ID grids has to be sufficiently far from every planet, 
so that the flux of angular momentum carried by the wave at the 
grid interface and the azimuthal dependence of the disk surface 
density are both small enough. An obvious test of our code is 
to check the dependence of the results on the positions of the 
interfaces between the two grids. 

We ran a bench of simulations like the one described in !3.4l 
with different inner radii for the 2D grid, but keeping constant 
the grid resolutions. Figure [TJshows the density profiles of the 
inner disk obtained at the end of the computations (16000 time 
units ~ 2500 orbits). It appears that the differences among the 
results are small, even in the cases where the planet is not very 
far from the interface (r p « 0.75). Only the cases with a transi- 
tion radius bigger than 0.5 show not negligible differences with 
the three other runs, which are remarkably similar to each other, 
in particular for what concerns the density profile of the inner 
edge of the gap. The case with interface at 0.5833 is clearly 
not accurate, but note that it is not unrealistic while the ID grid 
starts only at 4 Hill radii from the planet. 

All cases, however, give a quite consistent representation of 
the surface density profile of the inner disk, which is - on the 
contrary - very different from that obtained by using the 2D 
grid alone (compare with the dashed curve in Fig.|5}- 

This would not be the case if we had not implemented in 
the ID disk the exponential damping of the angular momen- 
tum carried by the density wave launched by the planet. For 
instance, Fig. [8] shows 3 of the simulations, recomputed by 
switching off the calculation of l|6}. One sees a kind of disconti- 
nuity at the interface, where the angular momentum deposition 
abruptly stops in the disk. This change implies a modification 



Fig. 8. Inner disk profiles at the end of simulations similar to 
those of Fig. but without the implementation of the wave 
damping algorithm A strong sensitivity on ^interface ap- 
pears. 



of the local equilibrium and of the shape of the density profile. 
Consequently, the results are strongly dependent on /interface- 
We note in passing that this also enlightens the importance for 
the gap structure of the flux of angular momentum carried away 
by the pressure supporte d wake. This flux is equivalent to the 
pressure torque studied in lCrida et al.l J2006). 

The outer interface has a much smaller influence on the 
disk's profile than the inner one, because it is further from the 
planet in the studied case. In Fig.|5] it appears that the outer 
interface corresponds to a change of sign in the second order 
derivative of the density profile. Fig.|5]shows that this is a real 
feature, specific to the chosen output time. Indeed, this figure 
shows the density profile at different times in the same simu- 
lation ; its evolution convincingly demonstrates that the second 
derivative of the density at the interface (/interface = 2.9167, 
right vertical dashed line) varies with time and can be non- 
zero. We also repeated the same simulation, moving the outer 
interface to /interface = 2.3167 (left vertical dashed line). The 
obtained surface density profiles overlap exactly with those of 
Fig- El at the corresponding time, so that they are not plotted. 
This test also demonstrates that the position of the outer inter- 
face has a negligible influence on the global disk evolution. 

The evolution of the disk has a strong influence on the 
type II migration of the planet. Figure [TO] shows the different 
migration rates in the simulations of Fig. [7] and in the simu- 
lation obtained with the 2D grid alone. Once again, the two 
first cases (/interface = 0.3333 and 0.4333) are remarkably sim- 
ilar, while the other ones show a slightly slower inward motion 
of the planet. However, the simulations done with the sole 2D 
grid, gives a migration that is sensibly faster, as expected due 
to the disappearance of the inner disk. 

In conclusion, the results presented in this section show that 
our method for coupling the 2D and ID calculations, not only 
ensures the conservation of the angular momentum, but also 
allows a robust (i.e. weakly dependent on the grid interface po- 
sition) modeling of the disk's structure and evolution. 
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Fig. 9. Outer disk profile close to the outer interface (/^interface = 
2.9167, marked by a vertical line) at 3 different times. The ver- 
tical line at r — 2.3167 marks the outer interface used in a sec- 
ond simulation, whose results are indistinguishable from those 
plotted in this figure. 
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Fig. 10. Type II migration of the Jupiter mass planet in the disk. 
The semi-major axis is plotted as a function of time for all the 
cases already discussed : various positions of the inner interface 
between the two grids, and no ID grid. 



4.2. Computational cost 

In this subsection, we show how the use of a ID grid allows the 
study of the whole disk for a negligible extra cost while the use 
of an extended 2D grid would be prohibitive. 

From a theoretical complexity point of view, the number 
of elementary operations for the computation of a time-step 
is proportional to the number of cells of the considered grid. 
Thus, it is obvious that the computation of the disk evolution 
in the ID grid is negligible with respect to the computation in 
the 2D grid. In addition, one has to consider the three opera- 
tions required to couple the grids, which also need to be done 
at every time-step. They are the filling of the ghost rings, the 
computation of SF/, and its spreading, and the change of T rS in 
the ID grid. They require computation on Nc 2D rings or less. 
As the size of the ghost area is usually negligible with respect 



to the size of the 2D grid, the computational cost of the cou- 
pling is also negligible with respect to the one of a time-step in 
the 2D grid. 

In the previous paragraph, we studied the number of ele- 
mentary operations to be computed during a time-step. In most 
codes, the length of the time-step is adapted to the conditions 
imposed by the required resolution and the existing perturba- 
tion. In general, it is de termined by the Coura nt-Friedrichs- 
Lewy (CFL) condition JStone & Normanl 1 19931 paragraph 
4.6). The limit most likely comes from the cells with smallest 
radii, where the angular velocity is the largest. Indeed, to avoid 
numerical instabilities, the time-step length 5t is set in order to 
avoid that information propagates more than one cell over one 
time-step ; thus one has 5t Q < 2n/N s , where Q is the angular 
velocity, which is about Keplerian in a gas disk. Consequently, 
denoting by R m \ n the inner boundary radius of the 2D grid, one 
has 6t oc R 3 ! 2 . This shows that extending the 2D grid toward 

min fc to 

the star shortens the time-step and slows down the simulation 
significantly. The use of the FARGO algorithm (Masset 2000) 
enables one to get rid of the mean angular velocity and to con- 
sider only the perturbed motion and the shear, but this leads 
to a similar conclusion, although the scaling of 5t with R m i n is 
generally different. 

The CFL condition in a ID grid is much less constraining. 
Consequently, the addition of a ID grid inside the inner edge of 
the 2D grid has no influence on the time step length, still deter- 
mined by what happens in the 2D grid. It is thus important to 
increase /interface as much as possible, to speed up the compu- 
tation. Some benchmarking confirmed this reasoning (see ap- 
pendix |E] for more detailed results). The accuracy loss of the 
computation with the increase of /^interface has been discussed in 
Subsect. f4.il so that an acceptable trade-off can be found. 



5. Astrophysical applications 

As mentioned in the introduction, the viscous evolution of the 
protoplanetary disk is thought to govern type II migration. Our 
code, which has been designed in order to correctly reproduce 
this evolution, while resolving the planet-disk interactions, is 
therefore very useful for studying any problem related to type II 
migration and the feedback exerted by the presence of planets 
onto the evolution of the disk. 

We present two problems for which this hybrid scheme is 
the tool of choice. 



5.1. Outward migrating planets 

A demonstration of how type II migration depen ds on the 
evolut ion of the disk has been provided by IVeras & Armi tage 
(2004). A disk evolving under its own viscosity accretes onto 
the central star, while spreading outward under the constraint of 
angular momentum conservation. Thus, at any instant in time 
there is a boundary in the disk, within which the radial motion 
of the gas in negative, and beyond which it is posi t ive - see 
the right panel of Fig. rTTland lLvnden-Bell & Pringld Jl974 . If 
a giant planet is located beyond this boundary it should move 
outward, as its mi gration has to follow th e local evolution of 
the disk (see also lLin & Papaloizoulll986l) . IVeras & Arm itage 
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(2004) showed this using a ID model, where the disk evolved 
under its own viscosity, the torque exerted in the planet-disk 
interaction had the form Q, and the effect of waves carrying 
an angular momentum flux was neglected. 

Using our code, we perform more precise simulations of 
this process. The interaction between the planet and the disk is 
simulated in 2D, and the effect of waves is taken into account. 
This implies that the ga p opened by the planet in the disk is less 
wide and deep tha n in lVeras & Armitaga (2004) representation 
(Crida et al. 2006). Moreover, the planet also feels a corotation 
torque, which is otherwise neglected in the analytic expression 
ofEq. 0. 

Inde ed we find that the situation is not as simple as illus- 
trated bv lVeras & Armitage] ( l2004]) . For instance, the left panel 
of Fig.^2 snows the result of a simulation in which the planet is 
initially put at r p = 1 .4, which is deeply in the outward spread- 
ing zone as shown by the right panel. Nevertheless, the planet 
migr ates inward, in what seems to be a runaway type III migra- 
tion (Masset & Papaloizou 2003). This is because the spread- 
ing disk forces enough gas to pass through the coorbital region 
of the planet s o that the planet feels a strong corotation torque 
( Mass e J200ll) and decouples from the disk evolution. 

We do obtain outward planet migrations, but only in spe- 
cific cases. For instance when we place the planet initially 
outside of the spreading disk (which is unrealistic), or if we 
hold it for a sufficiently large number of orbits, so that it can 
open a deep gap that effectively truncates the disk at its inner 
edge (which, in practice, places again the planet outside of the 
spreading disk). 

A detailed description of this mechanism and the com- 
prehensive exploration of the parameter space are beyond the 
scope of this paper. However, we think that this example is 
significant, because it shows that in reality the evolution of a 
Jovian mass planet is not an ideal type II migration. Because 
the gap is not extremely clean, the planet is not fully locked in 
the disk's evolution. The planet feels the global motion of the 
disk (accretion or spreading) but at the same time it feels also a 
non-negligible corotation torque. Only a code like the one that 
we have developed in this paper can simulate the two effects 
correctly and hence allows a quantitative study of the planet's 
evolution. 

5.2. Cavity opening 

If the migration of a giant planet is strongly influenced by the 
evolution of the disk, the presence of the planet also influences 
the evolution of the disk, with a non negligible feedback on its 
own migration. 

Perhaps the best example is that of the formation of an in- 
ner cavity in the disk. Once a planet has opened a wide and 
deep gap, the accretion time of the inner disk onto the star is 
smaller than the time scale of planet's migrati on, because o f 
the negative torque exerted by the planet ( Varniere et al. 2005). 
This leads to the disappearance of the inner disk. This in turn 
enhances the imbalance between the inner and outer torques 
felt by the planets, and accelerates the planet's migration. The 
presence of cavities are deduced in some protoplanetary disks, 
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Fig. 11. Left panel : migration of a Jupiter mass planet initially 
put on a circular orbit at r p = 1 .4. Right panel : Radial velocity 
in the gas disk when the planet is introduced. 

because the spectral energy distribution (SED) presents a lack 
of emissivity at the wavelengths correspondin g to hot gas i n 
the vicinity of th e star (see for instance ICalvetetalj l2005). 
iRice et al.1 d2003h suggest the SED of GM Aur could be ex- 
plained by a cavity maintained by a two Jupiter mass planet. 

Our algorithm clearly enables a more quantitative analysis 
of the cavity opening process than ever done before, for a low 
computational cost. As we have shown in Fig. [5] when the ID 
grid extends down to the real inner edge of the disk, the evolu- 
tion of the surface density of the inner disk is found to be slower 
than in a classical code with a truncated 2D grid, which makes 
the appearance of a deep cavity more difficult. Therefore, we 
think that the estimates of the planet mass required to achieve 
cavities of given depth, and the estimates of the lifetimes of 
these cavities, need to be revised. A detailed study of the cavity 
opening mechanism as a function of planetary mass, disk vis- 
cosity and aspect ratio is in progress and will be addressed in a 
forthcoming article. 

6. Conclusion and perspectives 

It is well known that the migration of giant planets (i.e. planets 
massive enough to open significant gaps in the disk's density 
distribution) is governed - or at least strongly affected - by the 
global evolution of the disk under its own viscosity. 

Usual simulation algorithms solve the hydrodynamical 
equations over a 2 dimensional polar grid that is truncated at 
an inner and outer radius to keep the computing time within 
reasonable limits. Thus, they cannot reproduce correctly the 
evolution of the disk, and no more the planet's migration. 

In this paper, we have shown how the use of a ID grid sur- 
rounding a classical 2D grid allows us to simulate the global 
evolution of the disk, while resolving the local planet-disk in- 
teractions with the accuracy of the usual algorithms. Coupling 
the two grids via a system of ghost rings, with special attention 
paid to the conservation of the angular momentum, leads to a 
smooth evolution of the disk from its inner radius (the trunca- 
tion radius or the surface of the star) to its outer edge, several 
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hundreds of AU away. This increases by far the accuracy of the 
simulation with essentially no additional CPU cost. 

Consequently, this algorithm is a tool of choice to properly 
simulate type II migration, and related problems. More gener- 
ally, it also enables the study of the effect of the presence of 
planets on the disk's global evolution, which in turn affects the 
migration of the planets in a feedback effect. 

In our simulations, we observe that the disks slowly dis- 
appear through the open boundaries of the ID grid. However, 
disks are believed to disappear rapidly after only a few million 
years, under the action of photo-evaporation. This phenomenon 
could easily be introduced in our code via a simplified prescrip- 
tion consisting in removing a fraction of the gas in each cell, 
depending on time and location. This will open the possibility 
of simulating the evolution of planets in a globally evolving gas 
disk, over the disk's lifetime. 
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Appendix A: Open boundary condition 

The "open" boundary condition, that allows outflow from the 
grid but not inflow is implemented as follows at the inner edge 
of the ID grid. 

In the innermost ring (number 0), the radial velocity is al- 
ways and the density is the one of the neighboring ring. In 
the latter (ring number 1), the radial velocity is set equal to the 
one of the following ring (number 2) if and only if it corre- 
sponds to outflow from the grid (negative radial velocity), and 
to otherwise. This gives : 

- 2(0) is set to 1(1). 

- Whenever v r (2) > 0, v r (l) is set to 0. 

- Whenever v r (2) < 0, v r {\) is set to v r (2). 

So, the first three rings are used for this computation, which 
explains the artifact observed in Figs.|^|5jEl and[S] 

Appendix B: Benchmarking 

We ran 3 simulations of a Jupiter mass planet initially placed 
at r p = 1, evolving in a gas disk for 1000 time units on an 
Intel Xenon 2.66 GHz processor. The results, summed up in 
table lB.ll confirm our reasoning of Sect. 14.21 

First, we see that with or without the ID grid, the comput- 
ing cost is about the same. 

Second, it is obvious that the computing time is not directly 
proportional to the number of cells, and is strongly increased by 
the higher shear in the vicinity of the star, otherwise the third 
simulation would not have been more that twice as long as the 
first one. 



Table B.l. Benchmarking results 



parameters of the 2D grid 


ID grid 


CPU time 


extension 


N r 


N s 


extension 


of process 


0.35 - 3.0 


165 


320 


none 


6167 s. 


0.35 - 3.0 


165 


320 


0.1167-20 


6008 s. 


0.1167-5.0 


294 


320 


none 


99494 s. 



Third, the second simulation was slightly shorter than the 
first one because with the ID grid, the density at the boundaries 
of the 2D grid does not tend to zero, the profile is smoother and 
consequently the CFL condition is less constraining there. 
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